Ecological and environmental factors affecting the risk of tick-borne encephalitis in Europe, 2017 to 2021

Background Tick-borne encephalitis (TBE) is a disease which can lead to severe neurological symptoms, caused by the TBE virus (TBEV). The natural transmission cycle occurs in foci and involves ticks as vectors and several key hosts that act as reservoirs and amplifiers of the infection spread. Recently, the incidence of TBE in Europe has been rising in both endemic and new regions. Aim In this study we want to provide comprehensive understanding of the main ecological and environmental factors that affect TBE spread across Europe. Methods We searched available literature on covariates linked with the circulation of TBEV in Europe. We then assessed the best predictors for TBE incidence in 11 European countries by means of statistical regression, using data on human infections provided by the European Surveillance System (TESSy), averaged between 2017 and 2021. Results We retrieved data from 62 full-text articles and identified 31 different covariates associated with TBE occurrence. Finally, we selected eight variables from the best model, including factors linked to vegetation cover, climate, and the presence of tick hosts. Discussion The existing literature is heterogeneous, both in study design and covariate types. Here, we summarised and statistically validated the covariates affecting the variability of TBEV across Europe. The analysis of the factors enhancing disease emergence is a fundamental step towards the identification of potential hotspots of viral circulation. Hence, our results can support modelling efforts to estimate the risk of TBEV infections and help decision-makers implement surveillance and prevention campaigns.


Introduction
Tick-borne encephalitis (TBE) is a zoonotic disease which affects human and animal central nervous systems with mild to severe long-term sequelae, which may be fatal [1,2].It is caused by the TBE virus (TBEV), a Flavivirus with currently three main subtypes and two additional subtypes recently proposed [3].They circulate in nature among ticks, mostly those belonging to the Ixodes ricinus complex, and in several wildlife hosts.The three main subtypes circulating in the European Union and European Economic Area (EU/EEA) are the European (Eu), Siberian (Sib) and Far Eastern (FE) subtypes [3].The European subtype TBEV-Eu, the most common one, is mainly associated with the biphasic form of TBE which has no chronic forms and presents symptoms with severe neurological sequelae in up to 10% of human cases and a fatality rate ranging from 1% to 2% [1].Transmission to humans usually occurs after a tick bite, although food-borne infections after consumption of unpasteurised milk and dairy products from infected animals have been reported [4].
The geographical occurrence of TBEV is fragmented, with foci of infection (hotspots) that are difficult to identify and often vary in space and time [5].Reporting of TBE cases in the EU/EEA is compulsory in 19 countries, voluntary in four (Belgium, France, Luxembourg and the Netherlands) and 'not specified' in one country (Croatia) [6], with 2,000 to 4,000 total cases reported yearly [7].The European Centre for Disease Prevention and Control (ECDC) has reported increases in TBE incidence over the last years [6].Major changes have been observed not only in the total number of reported cases, but also in the spatial distribution of the areas with active virus circulation, with the emergence of new TBEV foci in previously non-endemic countries [8][9][10].
The disease is preventable by vaccination along with personal protective measures which reduce the contact with infected ticks.The implementation of prevention and mitigation actions by public health authorities requires detailed knowledge of the disease's distribution, which, in turn, needs a comprehensive understanding of the ecological factors driving the intensity of viral circulation and infection hazard.
In recent decades, growing attention has been devoted to assessing the factors driving TBEV circulation within the natural foci, with several studies aimed at identifying abiotic (e.g.[11][12][13]) or biotic (e.g.[5,[14][15][16][17]) covariates.These include the analysis of the complex interactions between vectors and key vertebrate hosts that are strongly connected to the features of their local habitats.It is often difficult to establish the precise ecological conditions that favour TBE emergence and spread, a challenge that is reflected in the wide range of covariates that have been reported in the existing literature.
Our aim was therefore to obtain a more accurate understanding of the relationships occurring between a set of ecological variables and the incidence of TBE in humans across Europe, and to select those with highest impacts.As the TBE range has recently expanded and may continue to spread westward [8][9][10], northward [18][19][20] and to higher altitudes [21], we performed our analysis at a continental scale, responding to the need for a broader and more comprehensive understanding of the ecological forces driving such changes.This knowledge provides essential input to modern modelling approaches, based on quantitative disease data and a set of relevant covariates, which aim to predict the spatiotemporal risk of disease occurrence and its potential future spread.

Literature screening: search strategy and selection criteria
We performed a comprehensive literature search on TBE covariates following the principles of extending the PRISMA approach to scoping reviews [22].Keywords were extracted from the MeSH database and EMBASE vocabulary, then integrated with text words found in relevant papers; see Supplementary Table S1 for the search strategy and keywords.The search was performed on 21 July 2020.We used the CAS STNext platform to search the MEDLINE, EMBASE, BIOSIS, SCISEARCH and CABA databases.We also searched SCOPUS (via Elsevier) by adapting the search strategy to the database-specific characteristics.We included in our review primary research studies (i.e.studies generating new data), modelling studies proposing quantitative analysis using explanatory variables and data collections with abstract and full-text documents available in English, published after 1 January 2000.We excluded studies with no data or with duplicated data (patents, editorials, letters, modelling studies with no data).We also excluded records with no denominator, no identified reference population, unavailable fulltexts, and those that referred to data older than 2000 or were gathered outside the European Union and European Economic Area (EU/EEA).See Supplementary Figure S1 for the PRISMA flow diagram.

Four collaborators (AR, GM, LB, VT) independently evaluated potentially relevant records based on titles
What did you want to address in this study?During the last decades, the number of tick-borne encephalitis (TBE) cases reported in Europe has increased, making TBE a growing concern for public health.It is difficult to identify TBE risk areas, as the circulation of the TBE virus depends on the interplay between numerous environmental and ecological conditions.Our aim was to summarise all the different aspects that enhance TBE spread and identify the main forces that affect the distribution of TBE human infections in Europe.

What have we learnt from this study?
TBE is a seasonal disease, dependent on tick abundance and activity.We found that TBE spread is favoured by the presence of key animal species, such as deer and rodents, in forested areas.We also discovered that specific climatic conditions, such as high precipitation during the driest months of the year, cold winters, small daily temperature variations and a steep decrease in late summer temperatures, increase the risk of TBE infections in humans.

What are the implications of your findings for public health?
The identification of all the environmental and ecological aspects that are influencing the risk of TBE across Europe is fundamental for the rapid assessment of potential TBE outbreaks.Hence, this study will be used to inform future risk mapping efforts in Europe and in the long run improve the targeting of prevention and control measures.and abstracts.We then retrieved and read full texts of selected articles to assess their eligibility according to our inclusion and exclusion criteria and screened the references of selected publications to check for further sources of literature.We also added relevant articles published after the literature search was performed, up to 31 December 2021, by carrying out an additional search on PubMed.Finally, we used a pre-piloted data extraction spreadsheet to create our literature-based dataset and we selected covariates adopted in at least two articles for further analysis.

Epidemiological data
We analysed TBE case-based data provided by the European Surveillance System (TESSy) and released by ECDC.Each record included the date of disease onset, the importation status and the most probable place of infection.Coded values for variables with geographical information followed the European nomenclature of territorial units for statistics (NUTS).When available, the probable place of infection was provided at the NUTS-3 level, corresponding to small regions for specific diagnosis, according to Regulation (EC) No 1059/2003 [23].
Of all the TBE cases recorded in TESSy, we included only the laboratory-confirmed cases reported from 1 January 2017 up to 31 December 2021, since most countries did not report the place of infection before 2017.Patients infected outside their country of residence or whose location of exposure was unknown or provided at low spatial resolution, were excluded.We included countries that reported at least 10 cases between 2017 and 2021 and notified the place of infection at NUTS-3 level for at least 75% of the cases.The countries selected according to these criteria were: Czechia, Finland, France, Germany, Hungary, Italy, Lithuania, Poland, Slovakia and Sweden.We also included data reported from Austria although at a lower spatial resolution (corresponding to NUTS-2 regions, i.e. basic regions for the application of regional policies) as the spatial extent of NUTS-2 units in Austria is comparable to the NUTS-3 regions of the other countries.For each region i we computed the average annual TBE incidence Y i , expressed as the number of cases per 100,000 inhabitants, over the period 2017 to 2021.The total population in each spatial unit was extracted using gridded population count datasets (100 m spatial resolution) provided by WorldPop [24].We selected thirty-one covariates adopted in at least two articles for further analysis (Table 1).

Covariate data
We collected raw data from various sources to compute the covariates identified through literature screening.
The type of covariates considered were grouped into three different main categories, such as climatic, environmental and vertebrate host-related variables.
We used satellite images acquired by the moderate resolution imaging spectroradiometer and supplied by the National Aeronautics and Space Administration (NASA) with a resolution of 5.6 km as a source of land surface temperature and vegetation status as provided by the enhanced difference vegetation index (EVI).We downloaded the following products from the NASA Land Processes Distributed Active Archive Center: MOD11C1 Daily Land Surface Temperature and Emissivity [25], MOD11C3 Monthly Land and Surface Temperature and Emissivity [26] and MOD13C2 Vegetation Indices 16-Day [27].We computed cumulative precipitation data from the European Centre for Medium-Range Weather Forecast's fifth generation of European ReAnalysis (ERA5)-Land dataset and derived monthly time series of spatially enhanced relative humidity for Europe at 30 arc seconds resolution from ERA5-Land data [28].We calculated bioclimatic predictors following the formulae stated in the World Climate database [29] and computed averaged autumnal cooling and spring warming rates from 2017 to 2021 by applying a linear regression to the average daily temperature against the Julian day in the period 1 August to 31 October and 1 February to 30 April, respectively [30].
We extracted proportions of land cover classes from the 2018 Corine Land Cover database, with a resolution of 0.25 km.We calculated the total length of forest roads from raster maps of road density (km road per km 2 ), which were derived from linear road features extracted from OpenStreetMap datasets with a 1 km resolution.We derived estimates of snow and ice cover percentages from the 1 km consensus land-cover product [31].
Mean elevation was taken from the 1-km Global Multiresolution Terrain Elevation dataset [32].
To account for the distribution of hosts across Europe, we used 1-km data about the probability of presence of selected critical reservoir species (Apodemus flavicollis, Myodes glareolus) and a single variable that describes the probability of presence of cervid species (Dama dama, Cervus elaphus, Capreolus capreolus) that have been reported to be the most important amplifier hosts for I. ricinus with respect to other ungulate species [33].These variables were originally produced using spatial modelling techniques based on random forest and boosted regression trees [34,35].
We computed each covariate by averaging the raw values for the same spatial level as the available incidence data (NUTS-3 or NUTS-2).We also averaged covariate time series over the whole study period.

Statistical analysis
Firstly, we performed single-variable analysis aimed at investigating the association between each covariate x and TBE incidence Y i .The response variable Y i was log-transformed before analysis to normalise the distribution [36], and we included a random effect on the reporting country to consider potential differences among national notification systems.We defined second-order linear mixed models, one for each explanatory variable x, of the form: Where a 0 , a 1 and a 2 are the model coefficients, c is the random effect on the reporting country, and x i indicates the explanatory variable (Table 1).For each variable x i we tested both linear (L) and quadratic (Q) models.Quadratic models (Q) were selected as better models than linear (L) ones when the quadratic term proved to be significant (p value < 0.05).All explanatory variables, except those spanning an interval of 0-1 (EVI, land cover percentages, presence of hosts, rates of autumnal cooling and spring warming) were centred around their mean to avoid collinearity between linear and squared terms.
Afterwards, we used multiple linear regression to select the explanatory variables with the highest predictive power for TBE incidence (multivariable analysis).We built a full model considering all covariates with a significant (p value < 0.05) coefficient a 1 in the models previously described.For this subset, quadratic terms (coefficient a 2 ) were also included if significant in the single-variable analysis.All selected variables were examined for multicollinearity by computing Pearson's r pairwise correlation coefficients and variance inflation factors [37].Among highly correlated variables, we kept the ones with lowest Akaike information criterion (AIC) in single-variable analysis.We then computed all possible submodels and ranked them according to their AIC score.We finally selected the best parsimonious model with lowest AIC among a set of candidates with approximately equal performances (ΔAIC < 2) [38].Model assumptions were verified by checking the model's residuals for any pattern or dependency [39].We obtained p values according to the Satterthwaite method [40].All analyses were carried out using R v.4.1.2[41] and packages dplyr [42], exactextractr [43], raster [44], lme4 [45], lmerTest [46] and MuMIn [47].

Literature screening
After applying our selection criteria, we retrieved relevant information from 62 full-text articles (see Supplementary Figure S1 and Table S2 for the PRISMA flow diagram and list of articles).Most studies focused on central-eastern countries, such as Germany (16 studies) and Czechia (18 studies) (Figure 1A).The types of covariates considered were predominantly related to climate (46 studies), environment (37 studies) and competent and incompetent hosts (22 studies).We also included articles considering vector-host related data (12 studies), but as we selected covariates that were adopted in at least two articles for further analysis, none of the parameters used in such studies met this criterion (Figure 1B).The methodological approaches ranged from local surveys to more complex large-scale spatial models aimed at TBE risk assessment (Figure 1C).

Epidemiological and statistical analysis
In the period 2017 to 2021 a total of 12,289 confirmed cases with known place of infection were reported to ECDC from 371 NUTS-3 and nine NUTS-2 European regions from the 11 countries included in the study.The 4-year mean incidence across the considered NUTS ranged between 0.04 and 45.66, with an average (of all mean values) of 3.74 per 100,000 inhabitants (Figure 2).
The single-variable analysis proved that the distribution of the mean log-transformed TBE incidence transmission in Europe was significantly related to almost all factors, except for the total length of forest roads (For_length) (Table 2).
After checking for pairwise correlations (see Supplementary Figure S2 for the correlation matrix of covariates), we kept 22 covariates for multivariable analysis (see Supplementary Table S3 for the list of candidate models).Finally, eight covariates were selected in the best parsimonious model (Table 3).The model was characterised by a reasonably good fit (marginal R 2 = 0.28, conditional R 2 = 0.66, AIC = 921.56),higher than any of the fits obtained in single-variable analysis.Visual inspection of residual plots did not reveal any obvious deviations from normality.To better grasp how each predictor selected in the best model was related to the distribution of human TBE incidence in Europe, we computed conditional predictions for the log-transformed TBE incidence (log(Y i )) (Figure 3); all variables were kept at their average value, except for the one shown in each specific graph.
Our results show that higher TBE incidence in humans was linked to higher percentages of forested area and high precipitation in the driest quarter.Higher rates of autumnal cooling, a steep decrease in late summer temperatures, colder winters and smaller variations in daily temperatures values were also related to higher TBE incidence.Critical hosts species appear to have different impacts: disease incidence increased with the probability of presence of A. flavicollis, while it decreased in areas characterised by the presence of M. glareolus.We found a parabolic relationship between human incidence and the probability of presence of cervids (C.capreolus, D. dama, C. elaphus) (Figure 3).

Discussion
Tick-borne encephalitis is an increasing concern for European public health.The risk of infection depends on the co-occurrence of a set of ecological factors that have not been completely identified yet.Our literature screening revealed substantial heterogeneity in the selected studies.This diversity depends on the different goals of the studies, mostly focused on local investigations of TBEV in ticks [21,[48][49][50][51][52][53] and hosts [54][55][56][57][58].
Broader modelling studies assessing the geographical distribution of the pathogen are rarer and usually based on climatic predictors [59][60][61][62].Overall, we identified 31 covariates, and the single-variable analysis proved how TBE incidence was significantly affected by almost all of them.This result is in accordance with the available literature and provides additional confirmation of previous published analyses.Assessing the drivers shaping disease distribution is a fundamental step needed to successfully model disease risk.Eight factors proved to be the most effective for explaining the distribution of TBE incidence in Europe.
Firstly, it is essential to consider the presence of competent and non-competent tick-feeding hosts and the features of their habitat.Competent TBEV reservoir hosts are mainly small rodents and insectivores that support both virus circulation and feeding ticks, while non-competent hosts act as amplifiers of the vector population.All host-related variables were selected as relevant predictors in the best parsimonious model, underlining the importance of considering the presence of critical species when modelling the risk of emergence of new TBE hotspots.This is usually accomplished using local game animal density data as a proxy for host density [16,17,[63][64][65][66].However, it is difficult to retrieve such standardised information at the European scale, and these assessments are generally focused on non-competent hosts.In this study we used the probability of presence of rodents and cervids and validated their impact on the distribution of TBE incidence.These datasets could, therefore, serve as good predictors in future studies aimed at assessing the risk of disease outbreaks in vast geographical areas.
Rodents such as A. flavicollis and M. glareolus, the most common and widespread species inhabiting, often sympatrically, forested areas, play a pivotal role in the enzootic cycle of TBEV.They are well known susceptible hosts capable of transmitting the virus to the feeding ticks both systemically, developing viraemia, and non-systemically, via co-feeding [67].Interestingly, we found a positive relationship between TBE incidence and A. flavicollis, but a negative one with M. glareolus.One possible explanation could be the fact that M. glareolus acquires resistance to tick infestation [68], therefore hampering the co-feeding mechanism which is allegedly the most efficient mechanism contributing to TBEV circulation.
In addition to rodents, ungulates also play a major role in TBEV epidemiology [15,17,51,69,70] as they are able to amplify tick abundance by acting as hosts to adult stages and by moving them over long distances.At the same time, as non-competent hosts, they can divert tick bites from competent hosts (dilution effect), causing a decrease of TBEV prevalence in ticks after a certain threshold density is reached [15,51].Our results confirm this statement, as TBE incidence is lower in the regions characterised by low probability of presence of cervid species, then increases to reach a peak and finally decreases again in areas where the probability of occurrence is at its maximum.
The total proportion of forested areas (such as broadleaved, coniferous and mixed forest) was also found to be a good predictor for TBE incidence, with a positive impact on disease occurrence in humans.Forest areas provide suitable habitat and resources for ungulates, rodents and ticks, thus promoting their encounter rate and boosting the risk of occurrence of human TBE cases [55,[71][72][73][74][75][76].Moreover, human activity and behaviour can act in synergy with ecological and environmental factors by increasing the chances of exposure to infected ticks, as people engaged in recreational or occupational activities in forests are at increased risk of tick encounters and bites [77,78].The time spent in mixed forest for recreational purposes (of ≥ 10 h/week) has been positively associated with an increased TBE risk, and so were other activities such as harvesting forest foods and being employed as a forester or nonspecialised worker [77].
Tick-borne encephalitis is a seasonal disease, dependent on tick abundance and activity, which in turn is strongly affected by climatic conditions.We tested several variables related to temperature, precipitation and relative humidity in the full model and found four factors as the best predictors, namely the rate of autumnal cooling, the mean temperature registered in winter, the mean diurnal temperature range, and the total precipitation of the driest quarter.
At the continental scale, areas characterised by rapid temperature drops in late summer and early autumn are generally affected by higher values of TBE incidence, while at the local scale, the impact of daily temperature variations on the prevalence of TBEV in ticks based on field data showed contrasting results.For example, there was no evidence of any effect of a rapid autumnal temperature decrease on the minimum infection rate of nymphs in the following spring in a TBE focus in Germany [79].Such results were obtained by computing the decadal mean daily maximum air temperature in spring and autumn.On the other hand, the autumnal cooling rate (computed as in Randolph et al. [30]) proved to be a crucial ecological driver for cofeeding transmission of TBEV and for the maintenance of a TBE hotspot in northern Italy [54,69].Autumnal cooling plays a key role in TBEV epidemiology [47] as a steep decrease in late summer temperatures induces a behavioural diapause that favours a synchronous larval and nymphal activity the following spring, an event that is generally considered one of the most critical factors in TBEV transmission [30,54,69].
We hypothesise that the positive correlation between high TBE incidence and low winter temperature could be biased by the high incidence of cases registered between 2017 and 2021 in countries that exhibit low temperatures in winter, such as Austria, Czechia, Finland, Lithuania and Sweden.From an ecological perspective, this result can be explained by the fact that cold winter temperatures induce diapause in ixodid ticks, sheltering them from unfavourable climatic conditions and supporting their overwintering survival [80].On the other hand, TBE incidence decreases in regions characterised by strong daily temperature variations.Such changes in temperature may affect questing behaviour of ticks and thus the probability of contact with hosts and their survival [11,50].The total precipitation of the driest quarter is another key indirect factor that can influence tick behaviour and survival.Hence, higher precipitation might lead to lower tick mortality and continued tick questing during the driest months of the year, but also ensure that ticks in shelters survive to later activity periods [81].

Conclusion
TBEV distribution is shaped by the interplay of multiple climatic, environmental and ecological factors that exert a crucial role in the life cycle of ticks and TBEV circulation.Through our approach, we provided insights into the combination of covariates that appear to be crucial in affecting TBEV occurrence, defined their main data sources and established their interrelation with human TBE incidence at a large scale, considering the countries that notified TBE cases to ECDC at the highest possible spatial detail.The early identification of potential health threats derived from TBEV circulation is fundamental to improve timely detection and awareness of infectious disease events at the earliest stage of their emergence.Hence, this study could inform future modelling efforts aimed at assessing TBE risk across Europe and support competent authorities in deploying One Health integrated actions in existing and new potential risk areas.

Disclaimer
The views and opinions of the authors expressed herein do not necessarily state or reflect those of ECDC.The accuracy of the authors' statistical analysis and the findings they report are not the responsibility of ECDC.ECDC is not responsible for conclusions or opinions drawn from the data provided.ECDC is not responsible for the correctness of the data and for data management, data merging and data collation after provision of the data.ECDC shall not be held liable for improper or incorrect use of the data.

Ethical statement
This study is based on a literature search and is supported by ECDC epidemiological data.Covariate data have been retrieved from public sources.Ethical approval was not needed.
This is an open-access article distributed under the terms of the Creative Commons Attribution (CC BY 4.0) Licence.You may share and adapt the material, but must give appropriate credit to the source, provide a link to the licence and indicate if changes were made.
Any supplementary material referenced in the article can be found in the online version.
This article is copyright of the authors or their affiliated institutions, 2023.

Figure 1
Figure 1 Main characteristics of included studies on factors affecting the risk of tick-borne encephalitis, EU/EEA, 2000-2021 (n = 62)

Figure 3
Figure 3 Best model conditional predictions of factors affecting the risk of tick-borne encephalitis, EU/EEA, 2000-2021

Table 2
Results of the single-variable analysis on factors affecting the risk of tick-borne encephalitis, ordered by AIC, EU/EEA, 2000-2021 Significant coefficients are presented in bold.See Table1for a description of the predictors.